package gps

import "math"

//两点距离公式
func GetDistance(lat1, lng1, lat2, lng2 float64) float64 {
	radius := 6371372.0
	rad := math.Pi / 180.0
	lat1 *= rad
	lat2 *= rad
	lng1 *= rad
	lng2 *= rad
	theta := lng2 - lng1
	dist := math.Acos(math.Sin(lat1)*math.Sin(lat2) + math.Cos(lat1)*math.Cos(lat2)*math.Cos(theta))
	return dist * radius
}

//经纬度坐标转换到平面坐标
func MillierConvertion(lat float64, lon float64) (float64, float64) {
	var L, H, W, temp, mill, x, y float64
	L = 6381372 * math.Pi * 2 //地球周长
	W = L                     // 平面展开后，x轴等于周长
	H = L / 2                 // y轴约等于周长一半
	mill = 2.3                // 米勒投影中的一个常数，范围大约在正负2.3之间
	temp = math.Pi
	x = lon * temp / 180                           // 将经度从度数转换为弧度
	y = lat * temp / 180                           // 将纬度从度数转换为弧度
	y = 1.25 * math.Log(math.Tan(0.25*temp+0.4*y)) // 米勒投影的转换
	// 弧度转为实际距离
	x = (W / 2) + (W/(2*math.Pi))*x
	y = (H / 2) - (H/(2*mill))*y
	return x, y
}

//平面坐标转换到经纬度
func MillierConvertion_1(x float64, y float64) (float64, float64) {
	var L, H, W, temp, mill, lat, lon float64
	L = 6381372 * math.Pi * 2 //地球周长
	W = L                     // 平面展开后，x轴等于周长
	H = L / 2                 // y轴约等于周长一半
	mill = 2.3                // 米勒投影中的一个常数，范围大约在正负2.3之间
	temp = math.Pi
	lon = (x - (W / 2)) / (W / (2 * math.Pi))
	lon = lon * 180 / temp
	lat = -(y - (H / 2)) / (H / (2 * mill))
	lat = (math.Atan(math.Exp(lat/1.25)) - 0.25*temp) / 0.4
	lat = lat * 180 / temp
	return lat, lon
}

//求垂直点
func GetFootOfPerpendicular(x1, y1, x2, y2, x3, y3 float64) (float64, float64) {
	dx := x1 - x2
	dy := y1 - y2
	if math.Abs(dx) < 0.00000001 && math.Abs(dy) < 0.00000001 {
		return x1, y1
	}

	u := (x3-x1)*(dx) + (y3-y1)*(dy)
	u = u / ((dx * dx) + (dy * dy))

	return x1 + u*dx, y1 + u*dy
}

// WGS84坐标系：即地球坐标系，国际上通用的坐标系。
// GCJ02坐标系：即火星坐标系，WGS84坐标系经加密后的坐标系。Google Maps，高德在用。
// BD09坐标系：即百度坐标系，GCJ02坐标系经加密后的坐标系。

const (
	X_PI   = math.Pi * 3000.0 / 180.0
	OFFSET = 0.00669342162296594323
	AXIS   = 6378245.0
)

//BD09toGCJ02 百度坐标系->火星坐标系
func BD09toGCJ02(lon, lat float64) (float64, float64) {
	x := lon - 0.0065
	y := lat - 0.006

	z := math.Sqrt(x*x+y*y) - 0.00002*math.Sin(y*X_PI)
	theta := math.Atan2(y, x) - 0.000003*math.Cos(x*X_PI)

	gLon := z * math.Cos(theta)
	gLat := z * math.Sin(theta)

	return gLon, gLat
}

//GCJ02toBD09 火星坐标系->百度坐标系
func GCJ02toBD09(lon, lat float64) (float64, float64) {
	z := math.Sqrt(lon*lon+lat*lat) + 0.00002*math.Sin(lat*X_PI)
	theta := math.Atan2(lat, lon) + 0.000003*math.Cos(lon*X_PI)

	bdLon := z*math.Cos(theta) + 0.0065
	bdLat := z*math.Sin(theta) + 0.006

	return bdLon, bdLat
}

//WGS84toGCJ02 WGS84坐标系->火星坐标系
func WGS84toGCJ02(lon, lat float64) (float64, float64) {
	if isOutOFChina(lon, lat) {
		return lon, lat
	}

	mgLon, mgLat := delta(lon, lat)

	return mgLon, mgLat
}

//GCJ02toWGS84 火星坐标系->WGS84坐标系
func GCJ02toWGS84(lon, lat float64) (float64, float64) {
	if isOutOFChina(lon, lat) {
		return lon, lat
	}

	mgLon, mgLat := delta(lon, lat)

	return lon*2 - mgLon, lat*2 - mgLat
}

//BD09toWGS84 百度坐标系->WGS84坐标系
func BD09toWGS84(lon, lat float64) (float64, float64) {
	lon, lat = BD09toGCJ02(lon, lat)
	return GCJ02toWGS84(lon, lat)
}

//WGS84toBD09 WGS84坐标系->百度坐标系
func WGS84toBD09(lon, lat float64) (float64, float64) {
	lon, lat = WGS84toGCJ02(lon, lat)
	return GCJ02toBD09(lon, lat)
}

func delta(lon, lat float64) (float64, float64) {
	dlat := transformlat(lon-105.0, lat-35.0)
	dlon := transformlng(lon-105.0, lat-35.0)

	radlat := lat / 180.0 * math.Pi
	magic := math.Sin(radlat)
	magic = 1 - OFFSET*magic*magic
	sqrtmagic := math.Sqrt(magic)

	dlat = (dlat * 180.0) / ((AXIS * (1 - OFFSET)) / (magic * sqrtmagic) * math.Pi)
	dlon = (dlon * 180.0) / (AXIS / sqrtmagic * math.Cos(radlat) * math.Pi)

	mgLat := lat + dlat
	mgLon := lon + dlon

	return mgLon, mgLat
}

func transformlat(lon, lat float64) float64 {
	var ret = -100.0 + 2.0*lon + 3.0*lat + 0.2*lat*lat + 0.1*lon*lat + 0.2*math.Sqrt(math.Abs(lon))
	ret += (20.0*math.Sin(6.0*lon*math.Pi) + 20.0*math.Sin(2.0*lon*math.Pi)) * 2.0 / 3.0
	ret += (20.0*math.Sin(lat*math.Pi) + 40.0*math.Sin(lat/3.0*math.Pi)) * 2.0 / 3.0
	ret += (160.0*math.Sin(lat/12.0*math.Pi) + 320*math.Sin(lat*math.Pi/30.0)) * 2.0 / 3.0
	return ret
}

func transformlng(lon, lat float64) float64 {
	var ret = 300.0 + lon + 2.0*lat + 0.1*lon*lon + 0.1*lon*lat + 0.1*math.Sqrt(math.Abs(lon))
	ret += (20.0*math.Sin(6.0*lon*math.Pi) + 20.0*math.Sin(2.0*lon*math.Pi)) * 2.0 / 3.0
	ret += (20.0*math.Sin(lon*math.Pi) + 40.0*math.Sin(lon/3.0*math.Pi)) * 2.0 / 3.0
	ret += (150.0*math.Sin(lon/12.0*math.Pi) + 300.0*math.Sin(lon/30.0*math.Pi)) * 2.0 / 3.0
	return ret
}

func isOutOFChina(lon, lat float64) bool {
	return !(lon > 73.66 && lon < 135.05 && lat > 3.86 && lat < 53.55)
}
